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ABSTRACT 

We combine the GALFORM semi-analytical model of galaxy formation, which predicts 
the star formation and merger histories of galaxies, the GRASIL spectro-photometric 
code, which calculates the spectral energy distributions (SEDs) of galaxies self- 
consistently including reprocessing of radiation by dust, and artificial neural networks 
(ANN) , to investigate the clustering properties of galaxies selected by their emission at 
submillimctrc wavelengths (SMGs). We use the Millennium Simulation to predict the 
spatial and angular distribution of SMGs. At redshift z — 2, we find that these galaxies 
are strongly clustered, with a comoving correlation length of tq = 5.6 ±0.9 /i _1 Mpc for 
galaxies with 850/im flux densities brighter than 5 mJy, in agreement with observa- 
tions. We predict that at higher redshifts these galaxies trace denser and increasingly 
rarer regions of the universe. We present the predicted dependence of the clustering 
on luminosity, sub-millimetre colour, halo and total stellar masses. Interestingly, we 
predict tight relations between correlation length and halo and stellar masses, inde- 
pendent of sub-mm luminosity. 

Key words: galaxies :high-redshift - galaxies-evolution - cosmology: large scale struc- 
ture - submillimetre - methods:N-body simulations 



!_H \ 1 INTRODUCTION 



The discovery of a population of high-redshift galaxies 
selected by their emission at submillimetre wavelengths 
(submillimetre galaxies; SMGs) has opened a new win- 
dow on star formation in the h i gh re d shift Universe (e.g . 
Smail et all 1 19971; iHughes et al.l fl997l; iBarger et al l 1 19981 : 



Chapman etafl |2000| ; iBlain et al.1 l2002t ). The commonly 



held belief is that the submillimetre flux from these galax- 
ies is powered by prodigious star form ation rates, which 
can reach up ~ 5 00 - 1000 M a yr" 1 (|lvison et al.1 l200d : 
ISmail et al.1 (20021 : IChapman et al.l 120051 ). The star forma- 
tion is so intense that a substantial fraction of the mass 
of the present day descendants of SMGs, bright ellipticals, 
is thought to have been put in place du ring this phase 
l|Borvs et al.ll2005l ; lMTchalowski et al.ll2010l h A further con- 
straint on this picture would come from an estimate of the 
mass of the dark matter haloes which host SMGs. Measure- 
men ts of the clustering of SMGs have so far proved challeng- 
ing llScott et al.ll2002l: IBlain et alj|2004 iBorvs et all 120051 : 
IScott et al.l 12009 : IWeifi et al.l 12003 ) . This situation has re- 
cently improved with the launch of the Herschel telescope 
and will continue to get better with the commissioning of 
the new SCUBA-2 camera on the James Clerk Maxwell Tele- 



scope. In this paper we present predictions for the clustering 
of SMGs using a galaxy formation model set in the frame- 
work of the cold dark matter (CDM) cosmology. 

The self-consistent modelling of SMGs presents a num- 
ber of challenges. The sub-mm flux from a galaxy depends, 
often quite strongly, on a number of galaxy properties and 
parameters of the dust model, such as the star formation 
rate, the choice of the stellar initial mass function (IMF), 
the dust extinction (which is driven by the optical depth of 
the galaxy, which in turn depends on the metallicity of the 
cold gas and the size of the galaxy), the nature and composi- 
tion of the dust grains and the thermal equilibrium temper- 
ature reached by the dust grains when heated by starlight. 
Granato et al. (2000) introduced a hybrid model which com- 
bined a calculation of the star formation histories of galax- 
ies from the GALFORM semi-analytical galaxy formation model 
(Cole et al. 2000) with the spectro-photometric code GRASIL 
(Silva et al. 1998), which includes radiative transfer through 
a two-phase dust medium and a self-consistent prediction of 
dust temperatures. Using this model, a self-consistent cal- 
culation of the dust emission from galaxies can be made 
l|Baugh et al.ll2005l : lLa"cev et al.ll200al2010al h 

Constructing a galaxy formation model which can re- 
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produce the observed number counts of SMGs is relatively 
straightforward. It is more challenging to go a step further 
and to match the number counts and the redshift distribu- 
tion of SMGs at the same time. The task becomes much 
more difficult if, at the same time as matching the prop- 
erties of high redshift galaxies, the model is also required 
to reproduce obse rvations of the local galaxy population. 
iBaugh et all (2005) argued that it is only possible to achieve 
both of these goals by changing the slope of the IMF in 
episodes of star formation triggered by galaxy mergers. By 
adopting a top heavy IMF in starbursts, and by making the 
implied changes to the yield of metals and the fraction of gas 
recycled from dying stars, Baugh et al. were able to repro- 
duce basic properties of SMGs. At the same time, the same 
model matches the observed luminosit y function of Lyman - 
break galaxies from z — 3 to z = 10 (|Lacev et al.ll2010tj ). 
as well as being in good agreement with many observations 
of local galaxies. Nevertheless this model remains controver- 
sial and leads to conclusions which challenge the commonly 
accepted wisdom about SMGs. For example, Gonzalez et al. 
(2010b) investigated the nature of SMGs in the Baugh et al. 
model and found that the SMG phase is not responsible for 
the formation of a significant amount of long lived stars. 
Here we present a further test of the model by presenting 
predictions for the clustering of SMGs. 

To date, there are suggestions that SMGs are strongly 
clustered, with a correlation length that is substantially 
larger than that expected for the dark matter at the typical 
redshift of SMGs (e.g. Blain et al. 2004). However, the clus- 
tering measurements are currently noisy as a result of the 
small volumes surveyed, with this scatter being exacerbated 
by the strong clustering of the SMGs (Scott et al. 2002; 
2006; Borys et al. 2004; Weiss et al 2009). There is only lim- 
ited agreement between estimates of the angular clustering 
of SMGs, and poor agreement in turn between these mea- 
surements and the clustering inferred in three dimensions. 
The first results from the Herschel mission demonstrate the 
challenge of measuring the clustering of SMGs. Cooray et al. 
(2010) reported a detection of angular clustering in a sample 
of galaxies selected at 250 /im from the Herschel HerMES 
survey, which they estimate to have a similar redshift dis- 
tribution to the "classical" SMGs selected at 850 /im, while 
Maddox et al. (2010), on the other hand, found no evidence 
for angular clustering for a galaxy sample selected in a simi- 
lar way from the Herschel ATLAS survey. However, this sit- 
uation is likely to improved rapidly as the Herschel surveys 
increase in size and are analysed in more detail. Also, the 
SCUBA-2 camera is currently being installed at the JCMT. 
The Cosmology Legacy survey using SCUBA-2 will produce 
a map of 35 deg 2 at 850 fim, substantially bigger than the 
SHADES Half-Degree Extragalactic Survey. 

There are currently few predictions for the clustering of 
SMGs. van Kampen et al. (2005) compiled predictions for 
the angular clustering of SMGs from several groups. These 
calculations were phenomenological and did not attempt to 
predict the sub-mm flux from galaxies. The models were con- 
strained by hand to match the SMG number counts. Here 
we make a direct prediction of which galaxies satisfy the se- 
lection criteria to appear in an SMG sample. Gas dynamic 
simulations are currently unable to provide meaningful pre- 
dictions as the box sizes used are too small to predict clus- 
tering robustly beyond a scale on the order of a megaparsec. 



Furthermore, in many cases these calculations stop at high 
redshift (again due to the small box size) and so cannot 
be tested against the local galaxy population. By using a 
semi-analytical approach, the computational resources can 
be devoted to modelling the evolution of the dark matter 
component, allowing us to use a representative cosmological 
volume. 

In this paper we use the galform+grasil model to 
populate the Millennium Simulation of th e evolution of 
struc ture in a cold dark matter universe (jSpringel et al.l 
2005). This simulation occupies a volume of 500 h *Mpc on 
a side and contains more than 20 million dark matter haloes 
at the present day. The CPU required by GRASIL makes it 
impractical to compute the spectral energy distribution for 
every galaxy directly in the Millennium Simulation. Instead, 
we apply a novel technique based on artifical neural net- 
works (ANN) which we intr oduced in Paper I to populate 
the simulation with galaxies ([Almeida et al.ll2010l ). 

The paper is organised as follows. In Section 2 we give a 
brief summary of the GALFORM+GRASIL model and explain 
how the artificial neural network is implemented. We show 
how well the model can predict the submillimetre luminosity 
of galaxies in Section 3. In Section 4 we present the predic- 
tions for the spatial and angular clustering of SMGs. The 
dependence of the clustering on selected galaxy properties 
is explored in Section 5. Finally, in Section 6, we present our 
conclusions. 



2 MODEL 

Here we give a brief summary of the semi-analytical galaxy 
formation model, GALFORMfS I2.1|l . the spectro-photometric 
model used to compute galaxy SEDs, GRASIl(§ 12. 2[) and 
the artificial neural network (ANN) technique used to pre- 
dict spectral properties for large samples of galaxies. Fur- 
ther details a n d tes ts of this approach can be found in 
lAlmeida et ail (|2010l L 



2.1 The galaxy formation model: GALF0RM 

In this paper we use the GALFORM galaxy formation model 
to follow the fate of baryons in a ACDM universe. The gen- 
eral methodology behi nd sem i -analy tical modelling is ex- 
plained in the review by [Baugh ( 2006), and a more advanced 
overvi ew of galaxy formation physi cs is given bv | Benson! 
|2010l ). galform was introduced bv lCole et ail (|200Ct) . De- 
scriptions of subsequ e nt extensions t o the mode l are given in 
iBenson et aD lj2003l ): IBaugh et all (|2005h and iBower et al] 
|2006h . 

A summary of the model used in this paper, t hat of 



Baugh et a 



gn et a .1 (120051) . can be fo und in lLacev et aD ()2008l . 
)al) and lAlmeida et all (^Old ). where the processes mod- 



l2010ah and 

elled are described and the parameters used to specify the 
model are listed. Two changes from the original Baugh et al. 
set up are made as a result of the implementation of the 
model in the Millennium Simulation of Springel et al. (2005) . 
Firstly, the cosmological parameters of the Millennium are 
different from those adopted in the Baugh et al. model. E| We 



1 The Millennium Simulation adopts a flat ACDM cosmology 



Clustering of sub-mm galaxies 3 



have found that by adjusting the baryon density parameter 
from the Baugh et al. value of Sib = 0.045 to S7b = 0.033 
to give the same baryon fraction, f2b/O m , as used in the 
original model, we obtain similar predictions for the galaxy 
luminosity function to those obtained in the original model. 
Secondly, we use the merger histories of the dark matter 
haloes extracted directly from the Millennium, constru cted 
using the prescription described in iHarker et alJ <|2006h . 

An important feature of the Baugh et al. model, partic- 
ularly for the properties of galaxies selected by their dust 
emission, is the form of the stellar initial mass function 
(IMF) adopted in different modes of star formation. Bursts 
of star formation, which in this model are triggered by cer- 
tain types of galaxy merger, are assumed to produce stars 
with a top-heavy IMF, where diV/dlnm oc m~ x and x = 0. 
Bursts are initiated by all major mergers (i.e. those in which 
the mass in cold gas and stars of the accreted galaxy account 
for 30% or more of the primary's mass) and by minor merg- 
ers in which the accreted satellite makes up at least 5% of the 
primary's mass and where the primary is gas rich (defined 
as 75% of the primary mass being in the form of cold gas; 
these figures are model parameters). Quiescent star forma- 
tion in galactic disks is assumed to produce stars acc ording 
to a solar neig hbourhood IMF, the kennicuttl (|l983h IMF, 
with x = 0.4 for m < M© and x — 1.5 for m > M®. 

The adoption of a top-heavy IMF in starbursts is the 
key to reproducing the observed number counts and redshift 
distribution of faint sub-mm galaxies l|Baugh et all [2005; 
ISwinbank et alj|200ct ). While this choice is controversial, a 
variety of observational evidence suggests that in some en- 
vironments the IMF may have a higher proportion of high- 
mass star s than in the sol ar neighbourhood IMF (see the 
review bv lElmegreenl koOO). Moreover, the semi-analytical 
model is ideally placed to investigate the consequences for 
other predicted properties of assuming a top-heavy IMF in 
bursts. A number of predictions have been found to be in 
better agreement with observations following the use of dif- 
ferent IMFs in the burst and quiescent modes of star for- 
mation, such as the metalli cities of intra-cluster gas and of 
stars in early- type galaxies fNaga shima et al . 2005a b). The 
precise form of the IMF is not important so long as a higher 
proportion of high mass stars are produced than would be 
the case with a solar neighbourhood IMF. Similar predic- 
tions would be obtained for an IMF with a standard slope 
which is truncated below a few solar masses. With a larger 
fraction of massive stars produced relative to the Kennicutt 
IMF, more energy is radiated in the UV and larger amounts 
of dust are produced due to the enhanced yield of metals. 

In the next subsection we describe the GRASIL spectro- 
photometry code which generates a spectral energy distribu- 
tion for each galaxy across a wide range of wavelengths. GAL- 
FORM itself makes an independent calculation of the spectral 
energy distribution (SED) of starlight, including a mode l 
for dust extinction which is described in lCole et alJ l)2000t h 
This calculation gives similar results to those obtained with 



with a present-day matter density f! m = 0.25, a cosmological 
constant of Q/± = 0.75, a Hubble constant of h = _ffo/(100 km 
s -1 Mpc — 1 ) = 0.73 and a perturbation amplitude given by the 
linear rms fluctuation in spheres of radius 8h _1 Mpc of erg = 
0.9. The original Baugh et al. model also assumes a flat ACDM 
cosmology but with f! m = 0.3, S1a = 0.7, h = 0.7 and <rg = 0.93. 



GRASIL at optical wavelengths. The GALFORM calculation of 
the V-band luminosity weighted age and optical depth are 
used as inputs to the ANN. 

2.2 The spectro-photometric model: GRASIL 

To accurately predict the SEDs of galaxies, from the far- 
UV to the radio, we us e the spectro-photometric code 
GRASIL (|Silva et alJ 1 19981 ). This code computes the stel- 
lar emission, absorption and emission of radi ation by dust, 
and radio emission powered by massive stars (|Bressan et al.l 
2002). GRASIL carries out an accurate treatment of the ex- 
tinction and reprocessing of starlight by dust. 

The combinat i on of G ALFORM and GRASIL was described 
by iGranato et all ife oOO) and ha s been explo ited in a se- 
ries of papers (jBaueh et al.ll2005l : lLacev et al1l200Sl . l2010al ) . 
The semi-analytical model calculates the star formation and 
metal enrichment history for each galaxy, including the con- 
tribution from starbursts. GALFORM also predicts the scale- 
lengths of the disk and bulge components of each galaxy, as 
descri bed i n Cole et al. (2000) and tested by I Almeida et all 
l|2007l ) and iGonzalez et al] ((2009), an d the cold ga s mass 
(as compared against observations by I Power et al] (|2010t ) 
and iKim et al] l|201Ch ). The dust is modelled as a two- 
phase medium, with a diffuse component and dense molec- 
ular clouds. The mass split between these components is a 
model parameter. In the Baugh et al. model, 25 per cent of 
the dust is assumed to be in the form of dense clouds. Stars 
form within molecular clouds and esc ape on a time s cale t esc , 
which is another model parameter; in lBaugh et al.l l|2005l ). a 
value of t csc — 1 Myr is adopted in both quiescent and burst 
modes of star formation. The extinction of starlight by dust 
clouds depends on the star's age relative to the escape time. 
High mass stars, which typically dominate the emission in 
the UV, spend a significant fraction of their comparatively 
short lifetimes within molecular clouds. GRASIL calculates 
the radiative transfer of starlight through the dust and self- 
consistently solves for the temperature distribution of the 
dust grains at each point in the galaxy, based on the local 
radiation field. The temperature distribution of the grains is 
then used to calculate the dust emission. The composition 
and size of the dust grains are chosen to match the properties 
of the local ISM: a mixture of graphite and silicate grains, as 
well as polycyclic aromatic hydrocarbon (PAH) molecules. 
The effects of temperature fluctuations in very small grains 
and PAH molecules are taken into account. E mission from 
PAHs is calculated using the cross-sections of iLi fc Draind 
(2001). Radio emission from ionis ed HII regions and s yn- 
chrotron radiation is included as in lBressan" et al l l|2002l ). 

The GRASIL model has been calibrated against lo- 
cal observational data for normal and starburst galaxies 
llBressan et al.ll2002l : I Vega et al.ll2005l : IPanuzzo et~al]|2007l : 
ISchurer et al] 12009). A limitation of GRASIL is that it as- 
sumes axisymmetric distributions for the gas and dust in 
starburst galaxies. There is observational evidence of more 
complex geometries and extraplanar dust in some galax- 
ies (see for example Engclbracht et al. 2006). This could be 
problematic if this dust absorbed and emitted a significant 
fraction of radiation. However, there is little observational 
evidence for this. Furthermore, observations of nearby star- 
bursts reveal that most of the absorption and emission of 
radiation by dust takes place in a compact region of size ~ 1 
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kpc or less. GRASIL has been shown to accurately predict the 
SEDs of both quiescent a nd starburst galaxies IjSilva et al.l 
Il998l : iBressan et alll2002h . 



2.3 The artificial neural network approach to 
predicting galaxy luminosities 

The GRASIL code provides an accurate calculation of the 
absorption and reemission of radiation by dust, predicting 
the SED of a galaxy from the far-UV to radio. However, 
GRASIL is extremely CPU intensive, requiring several min- 
utes to compute the SED for a single galaxy, which pro- 
hibits its d irect application to e xtremely large numbers of 
galaxies. In I Almeida et all (|201Cf ) we introduced a new tech- 
nique based on artificial neural networks (ANN) which can 
be used to rapidly predict SEDs using a small set of galaxy 
properties as input, once the ANN has been trained on a rel- 
atively small number (~ 2000) of galaxies with SEDs com- 
puted using GRASIL. We demonstrated that, in the majority 
of cases, this method can predict the luminosities of galaxies 
to within 10 per cent of the values computed directly using 
GRASIL. We employ the same approach in this paper. The 
general methodology behind the ANN is set out in detail in 
Almeida etafl feoich . so we give only a brief summary here. 
Silva et al.l Q2010j) recently published a complementary ap- 
proach in which the explicit calculation of emission by dust 
within GRASIL is replaced by an ANN. 

Artificial neural networks are mathematical models de- 
signed to replicate the behaviour of the human brain. They 
are similar to their biological counterparts in the sense that 
ANNs consist of simple computational units, neurons, which 
are interconnected in a network. The neurons are usually or- 
ganized in layers: an input layer, one or more hidden layers 
and an output layer. Each neuron has a weight associated 
with it. In this paper we will use a multilayer, feed- forward 
network. The ANN is trained using a sample of galaxies 
for which GRASIL has been run to compute spectra. Dur- 
ing the training process, the neural network is presented 
with a set of inputs, comprised of selected galaxy proper- 
ties, and associated outputs, in our case the luminosity at 
different wavelengths. The network weights are adjusted in 
order to reproduce, as closely as possible, the desired out- 
put from the given set of inputs. In summary: (i) we start 
with an untrained net (random weights); (ii) determine the 
output for a given input; (iii) compute the discrepancy or 
error between the predicted and the target output; and (iv) 
adjust the weights in order to reduce this error. To adjust 
the weight s we use the resilient bac kpropagation learning 
algorithm jRiedmiller fc Braunlll993l l. 

As in I Almeida et all l|201Ch . we use 12 galaxy proper- 
ties predicted by GALFORM as input to the ANN: the total 
stellar mass, the stellar metallicity, the unextincted stellar 
bolometric luminosity, the disk and bulge half-mass radii 
and the circular velocities measured at these two radii, the 
V-band weighted age, the optical depth of dust extinction in 
the V-band, the metallicity of the cold gas, the mass of stars 
formed in the last burst and the time since t he sta rt of the 
last burst of star formation. lAlmeida et al. (2010) showed 
that the performance of the ANN is greatly improved if we 
predict only one output property, the luminosity in a single 
band-pass, instead of predicting the full SED of the galaxy 
(which typically covers 456 wavelength bins for a standard 



32 



26 



24 



850 /im 




Figure 1. Comparison between the ANN predicted luminos- 
ity, ^predicted! and the true luminosity calculated using GRASIL, 
itruo, at 850 /im in the observer frame for galaxies at z = 2. 
Note that a flux of 1 mjy at z = 2 corresponds to a luminosity 
L„ = 4.75 X 10 31 /i _2 ergs _1 Hz _1 . The shading shows the distri- 
bution of galaxies in the training sample. In the inset, we plot 
the error distribution of the predicted luminosities, as given by 



^predicted 



/Lt 



1, normalized to unit area. 



GRASIL SED). We then train the ANN separately for each 
band required. Here, we follow the same approach: we train 
two separate networks, one for each of the sub-mm wave- 
lengths at which we want to predict luminosities (450 /im 
and 850 /im). The network configuration adopted has 12 neu- 
rons in the input layer, two hidden layers with 30 neurons 
each and one output neuron. We also found, in order to 
maximise the accuracy of the ANN predictions, that it was 
necessary to train the ANN at each redshift of interest, and 
to train separately for galaxies whose star formation is dom- 
inated by starbursts or quiescent star formation in disks. 



2.4 The performance of the ANN at sub-mm 
wavelengths 

We now demonstrate the how well the ANN performs when 
predicting galaxy luminosities in the sub-mm . Test at other 
wavelengths were presented in lAlmeida et ail (|201Ch . 

In Fig. [T]we plot the comparison between the observer 
frame luminosity in the SCUBA 850 ^m band predicted by 
the ANN for z = 2 galaxies and the true values calculated 
directly using GRASIL. In this plot we include all galaxies 
regardless of their classification as quiescent or starburst 
galaxies. Fig. Q] shows that there is excellent agreement be- 
tween the luminosities predicted by the ANN and the true 
values, with most of the predicted luminosities being within 
10 per cent of the GRASIL result (inset) . Some statistics quan- 
tifying the error distribution at different redshifts are sum- 
marized in Table [T] for galaxies brighter than 1 mjy in the 
corresponding band (either 450 or 850 pirn). Here, the root 
mean squared logarithmic error, el, is defined by: 
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450 fim 




850 fim 


Redshift 


Sample 


el 


^ e | < 10 per cent 


el 


^ e < 10 per cent 


z = 0.1 


Quiescent 
Burst 


0.04 
0.07 


97.1 
88.3 


0.04 
0.13 


98.6 
87.0 


z = 0.5 


Quiescent 
Burst 


0.04 
0.05 


96.9 
94.8 


0.04 
0.05 


98.5 
96.7 


z = 1 


Quiescent 


0.06 


94.9 


0.04 


97.0 


Burst 


0.08 


90.1 


0.07 


93.2 


z = 2 


Quiescent 


0.04 


97.2 


0.04 


97.5 




Burst 


0.05 


92.1 


0.05 


95.1 


z = 3 


Quiescent 


0.05 


96.0 


0.03 


97.8 




Burst 


0.05 


93.4 


0.03 


98.2 


z = 4 


Quiescent 


0.04 


97.0 


0.03 


98.5 




Burst 


0.05 


93.2 


0.03 


97.9 



Table 1. Statistics of the error distribution associated with the ANN prediction of 450 and 850 fim observer- frame luminosities at 
selected redshifts. The statistics are computed using only galaxies with sub-mm fluxes brighter than 1 mjy in that band. Column 1 gives 
the redshift, column 2 specifies whether the galaxy sample is made up of galaxies forming stars quiescently or starbursts; columns 3 and 

4 give el (the root mean squared error denned by Eq. [TJ and -P| e |<io per cent (percentage of galaxies with predicted luminosities within 
10% of the true value), respectively, for the 450 fim predictions. For 850 fim selected galaxies, the same information is shown in columns 

5 and 6. 



£l = y l/n^[ln(Lp rcdictcd /L truc )] 2 , (1) 

where n is the number of galaxies considered. The quan- 
tity -P| e |<ioperccnt is defined as the percentage of galaxies 
with predicted luminosities which lie within 10 per cent of 
the true values. For quiescent SMGs, we are able to repro- 
duce the luminosities of more than 95 per cent of galaxies 
with an accuracy of 10 per cent or better, for the redshift 
range consider e d, at both 450 fim and 850 fim. As shown by 
lAlmeida et al.l (2010) , the performance of the ANN for burst 
galaxies is somewhat poorer, which is a consequence of the 
wide range of spectra seen in bursts and the difficulty the 
ANN experiences in reproducing this variety. Nonetheless 
the technique returns more than 90% of predicted sub-mm 
luminosities within 10 per cent of the true values. It should 
be noted that at z = 4, the highest redshift considered, 
the observer-frame 850 fim luminosity probes the rest-frame 
170 fim, which is approaching the peak in the dust emission 
spectrum (typically around 100 /im). 

One important feature of the error distribution is shown 
in the inset of Fig. [T] The distribution of luminosity errors 
predicted by the ANN appears to be Gaussian. Furthermore, 
we find that there is no correlation of the error with luminos- 
ity or other galaxy properties. This suggests that any sample 
of SMGs built using the ANN method will have errors which 
are decoupled from the structural and photometric proper- 
ties of the galaxy sample. 



3 CLUSTERING OF SUB-MM GALAXIES 

The clustering of galaxies is an important constraint on the 
masses of their host dark matter halos, and hence on theo- 
retical models of galaxy formation, as it depends upon how 
various physical processes vary with halo mass. In this sec- 
tion we present the model predictions for the clustering of 
galaxies selected by their flux at sub-mm wavelengths. In 
§3.1, we contrast the spatial distribution of galaxies with 
that of dark matter haloes. In §3.2, we define the two-point 



spatial correlation function. We demonstrate that our clus- 
tering predictions are insensitive to errors in the galaxy lu- 
minosities predicted by the ANN in §3.3. We present the 
predictions for clustering in real space and redshift space 
in §3.4 and §3.5 respectively. The evolution of the correla- 
tion function is presented in §3.6 and the angular correlation 
function is shown in §3.7. 

3.1 The spatial distribution of SMGs 

Before discussing the predictions for the two-point correla- 
tion function of SMGs, we first gain a visual impression of 
their spatial distribution. Fig. [2] shows dark matter haloes 
and sub-mm galaxies in a slice taken from the Millennium 
Simulation. The slice measures 100 h' 1 Mpc across and 
20 k' 1 Mpc thick in comoving units. The upper panels show 
haloes and galaxies at z = 0.5 and the lower panels show 
them at z = 2. Dark matter haloes are shown by the blue 
shading. The intensity of the shading is proportional to the 
total halo mass within each pixel. We show sub-mm galaxies 
selected at 450 fim and 850 fim with fluxes brighter than 1 
mjy and 5 mjy. At a given flux limit, the 450 fim sources 
are more numerous than the 850 fim sources. The 850 fim 
sources brighter than lmjy tend to trace out the more mas- 
sive dark matter haloes and hence are expected to be biased 
tracers of the dark matter distribution. The full Millennium 
box (500ft _1 Mpc across) subtends an angle of 7.5 degrees at 
z — 2. To put this into context, we note that the SCUBA- 
2 Cosmology Legacy Survey (SCLSJ3 aims to map around 
35 deg 2 at 850 fim in patches up to 3 deg across, and 1.3 deg 2 
at 450 fim in regions up to 0.5 deg across. The nominal 5a 
flux limits will be 3.5 mjy at 850 fim and 2.5 mjy at 450 fim. 
However, source confusion may result in the flux limits for 
reliable source identification being somewhat brighter than 
this; using the standard 20 beam s per source criterion for 
confusion (e.g. lLacev et al.ll2008l ) together with the counts 

2 http:/ /www. jach.hawaii.edu/JCMT/surveys/Cosmology.html 
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Figure 2. Image of the simulated spatial distribution of SMGs and dark matter haloes at z = 0.5 (top panels) and 2 = 2 (bottom panels). 
The panels display a slice 100 /i -1 Mpc wide with a depth of 20 h~ 1 Mpc (in comoving coordinates). The width of this slice corresponds 
to angular scales of 4.3° and 1.6° at z = 0.5 and 2 = 2 respectively. The dark matter haloes are plotted in blue, with the darker shading 
corresponding to regions of higher projected halo mass density. Galaxies with S^gOpm and Ssso^m ^ 1 m Jy are represented by the orange 
dots. Brighter galaxies with S ^ 5 mJy are shown by the black dots. The left hand panels show galaxies selected by their 450 fim flux, 
while the right hand side show galaxies selected by their emission at 850 fim. 



predicted by the model, we expect confusion to become im- 
portant around 3.0 mJy at 850 /im and 5.3 mJy at 450 /im. 



3.2 The two-point correlation function 

To quantify the clustering of the galaxy distribution we use 
the two-point correlation function, £(r), which gives the ex- 
cess probability, compared with a random distribution, of 
finding two galaxies at a separation r: 

5P{r) = n 2 [l + Z(r)]5V 1 SV 2 , (2) 

where n is the mean space density of galaxies and the 8V\ 
are elements of volume. If f (r) > 0, then galaxies are more 



clustered than a random distribution. On the contrary, if 
galaxies have a tendency to avoid one other, then £(r) < 0. 

The two-point correlation function of galaxies is shaped 
by two main factors, which play different roles on different 
scales. On large scales, the form of the correlation function 
is controlled by the clustering of galaxies in distinct dark 
matter haloes (referred to as the two-halo term), and the 
galaxy and dark matter correlation functions have similar 
shapes but differ in amplitude (for an illustration of this see 
Angulo et al. 2008). On smaller scales, up to the size of the 
typical haloes which host galaxies, the form of the correla- 
tion function is driven by the number and radial distribu- 
tion of galaxies wit hin the same dark matter halo, (called 
the one-halo term) (jBenson et aUl200d : ISeliaklbOOOT ). 
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Figure 3. The impact of errors in the ANN predicted 850/^m 
luminosities on the form of the real-space two-point correlation 
function. The correlation function is measured for galaxies se- 
lected with fluxes brighter than 1 mjy at z = 2 (£ prc d)- The 
galaxy fluxes are then perturbed using the error distribution of 
the ANN 850 fim luminosities and the correlation function is re- 
measured for the new sample of galaxies brighter than 1 mjy 
(£pcrt)- The plot shows the maximum deviation of f pe rt from 
£prcd on constructing 20 different perturbed samples, expressed as 
a percentage. The spike around log(r/h — 1 Mpc) f» —0.8 is caused 
by noise due to the small number of galaxy pairs at that separa- 
tion. 



We calculate the two-point correlation function of sub- 
mm selected galaxies, using both real and redshift space 
coordinates. We measure the correlation function using the 
standard estimator (e.g. Peebles 1980): 



£(r) = 



(DD(r)) 



±7V gal nAF(r) 



- 1, 



(3) 



where (DD(r)} is the number of distinct galaxy pairs with 
separations between r and r + Ar, jV ga i is the total number 
of galaxies, and n is the mean number density of galaxies. 
A V(r) is the volume of a spherical shell of radius r and thick- 
ness Ar. We are able to compute the volume of this shell 
analytically since we are dealing with galaxy pairs within a 
periodic simulation volume. The clustering signal is gener- 
ated in redshift space using the distant observer approxima- 
tion, by electing one axis to be the line of sight direction, 
and adding the suitably scaled peculiar velocity of the galaxy 
along this axis to its comoving position. 

3.3 The impact of luminosity errors on the 
predicted correlation function 

We now look at the impact of errors on the ANN-predicted 
luminosities on the amplitude and shape of the two- 
point correlation function. First, the correlation function, 
£ P red(?"), is measured using the ANN-predicted observer- 
frame 850 /jm luminosities for a sample brighter than some 
flux limit, in this case 1 mjy. Next, these luminosities are 



perturbed by drawing an offset from the distribution of er- 
rors expected for this band (see inset of Fig. [T] and Table [TJ. 
A new flux limited sample is constructed, which will contain 
some galaxies which were not included in the initial, un- 
perturbed sample, because their fluxes have been boosted. 
Moreover some galaxies which made it into the original sam- 
ple will no longer be included after their luminosities have 
been perturbed. We then repeat the measurement of the 
correlation function for this new sample, resulting in the es- 
timate £ P ort(7"). Comparing ^pcrtM with £ pre d(r) provides an 
estimate of how £(r) is affected by the errors in the ANN- 
predicted luminosities. In Fig. Owe plot the maximum devi- 
ation of the ratio £ pe rt('")/£pred('") using 20 different £ per t(f) 
measurements for galaxies at z — 2 (i.e. after perturbing the 
galaxy luminosities). This plot shows that £ pal t{r) / £ ple d(r) 
differs from unity by at most 4%, indicating that the clus- 
tering predictions are essentially unaffected by the errors in 
the ANN luminosities. At z = 0.1, where the performance 
of the ANN is poorer for galaxies undergoing a burst of star 
formation than it is at z — 2, the ratio £ pe rt(r) / £pred{r) still 
deviates from unity by less than 10%. 

3.4 The real space correlation function 

In real space, the cartesian coordinates of the SMGs within 
the simulation box are used to calculate pair separations. 
Fig. U (top panel) shows the real space correlation function 
for sub-mm selected galaxies, for both A = 850 /im and A = 
450 /jm, at redshift z — 2. The transition between the one- 
halo term and two-halo term occurs around r ~ 0.6/i _1 Mpc. 
In this plot we consider samples of galaxies selected to be 
brighter than 1 mjy or 5 mJy at both wavelengths. The 
black line shows the correlation function of the dark matter, 
£dm, which was measured using a randomly chosen subset 
of 10 6 dark matter particles out of the ten billion particles 
in the Millennium Simulation. 

As can be seen in Fig. [4] galaxies do not trace the mass 
distribution in the Universe, because t he efficiency of g alaxy 
formation depends on halo mass fe.g. lEke et alj |2004). The 
difference between the clustering of galaxies and the un- 
derlying dark matter is quantified in terms of the galaxy 
clustering bias, b: 

& w=(r4r) 1/2 - < 4 > 

Udm [r)J 

Numerical simulations have demonstrated that the galaxy 
bias is a function of scale (Smith, Scoccimarro & Sheth 2007; 
Angulo et al 2008) . This scale dependence weakens on large 
scales, and the galaxy bias is typically approximated as a 
constant. We plot the bias of SMGs in the lower panel of 
Fig. 2] The plot shows that at z = 2 the bias factor is gen- 
erally greater than unity, approaching a roughly constant 
value of b ~ 1.8 for r > 2ft _1 Mpc for all of the samples 
shown. In the case of galaxies selected at 450 /im, there is 
a small but clear difference in the bias predicted for bright 
and faint samples, with the bright galaxies being the more 
strongly clustered. At 850 /jm the distinction is less clear, 
due in part to the relatively low number density of galaxies 
in the bright sample, which results in a noiser prediction. 

The effective bias parameter on large sca le s can 
also be estimated analy tically (|Mo fc White! 1 19961 ; 
ISheth. Mo fc Tormenl |200lh , using the mass function 
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Figure 4. The real space two-point correlation function of sub- 
mm selected galaxies in comoving coordinates at z = 2. The blue 
and orange dashed lines show the predicted correlation functions 
for galaxies with fluxes brighter than 5 mjy at 450 and 850 fj,m 
respectively, while the green and red lines show the corresponding 
results for a flux limit of f mjy. The black line shows the correla- 
tion function of dark matter. The errorbars show the la Poisson 
errors derived from the number of pairs in each bin of radial sep- 
aration, and are only shown for those samples with the lowest 
total number of galaxies. In the bottom panel, we plot the galaxy 
bias, b(r), as a function of scale for the four samples of sub-mm 
selected galaxies. This is obtained by taking the square root of 
the ratio of the galaxy correlation function to the measured dark 
matter correlation function. 



of haloes which host sub-mm galaxies, N(z,M), (i.e. the 
product of the space density or mass function of dark 
matter haloes and the halo occupation distribution of 
SMGs) and the bias factor as a function of the halo mass 
b(z,M) (e.g. iBaugh et al.lll999l ): 



bea(z) = 



J M N(z, M') b(z, M') d In M' 
f M N(z,M>) din M> 



(5) 



The integrals are taken over the full range of halo masses, 
with N(z,M) = for haloes which do not host SMGs. 
To compute b(z, M) we us e the prescription outlined by 
ISheth. Mo fe Tormenl (|200ll ). For galaxies at z = 2 with 
•S^Oum ^ 5 m Jyi we nn d an effective bias of 6 c ff = 2.3, and 
b c s = 2.1 for sub-mm galaxies with Siso^m iZ 5 mjy. These 
values are slightly larger than those estimated from the sim- 
ulation using Eq. [4] In fact, similar differences between the 
analytical approach and simu lations have been observed in 
other studies (see for example lGao et al 1 l2005l ; lAngui o et al.l 
l2009h . 

On small scales, Fig. [4] shows that the effective bias 
takes on a range of values. The clustering on these scales is 
driven by the typical number of galaxy pairs within a com- 
mon halo. The faint sample of galaxies selected at 450/im 
displays the strongest clustering on small scales. This sample 
contains the largest number of pairs within common haloes. 
In the case of the bright sample at 850/^m, the low number 



density of galaxies makes it difficult to measure the correla- 
tion function on small scales. 

A convenient measure of the strength of clustering for 
different galaxy samples is provided by the correlation length 
ro, which we can define in a robust way as the pair separation 
at which the correlation function becomes unity: 



«n>) = l. 



(6) 



Applying this definition to model galaxies selected at 
850 /im, we find ro = 5.6 ± 0.9 /i~ 1 Mpc for Sssoum ^ 5 mjy, 
and ro = 5.38 ± 0.02/i _1 Mpc for a fainter sample with 
Sgso/jm ^ 1 mjy. We also find that 450 /im selected galax- 
ies are less clustered than 850 ^m selected galaxies at the 
same flux limit: we obtain ro = 5.38 ± 0.02/i _1 Mpc and 
4.99 ± 0.01 h~ 1 Mpc respectively for S^o^m ^ 5 mjy and 
1 mjy. SMGs with Sisofim ^ 5 mjy display a similar two- 
point correlation function to that of Sssoijm ^ 1 mjy (as we 
will see later, this is mainly a consequence of the fact that 
the median S450m™/>S'850m™ colour is approximately 3). 

An alternative way to define the clustering length ro, 
commonly used in observational studies, is to fit the corre- 
lation function with a power law: 



For optical galaxies at z ~ 0, this is found observationally 
to provide a good fit for 0.1 < r < 10/i~ 1 Mpc with 7 close 
to -1.8 (e.g. iNorberg et all l200ll : IZehavi et all 120051 ). The 
two definitions of ro (Eqs. [6] and are obviously equiv- 
alent only if £(r) really is a power law. If £(r) actually 
has a more complicated dependence on r, then the value 
of ro obtained by fitting a power law will depend on the 
range of r over which the fit is performed, and on the er- 
rors on the measurements at different r. If we fit £(r) for 
our model SMGs at z = 2 with a power law over the range 
0.1 < r < 10/i _1 Mpc, we find 7 = -1.8 ± 0.2 and a cor- 
relation length of ro = 6 rt 1 h~ Mpc for Sg^oiim 5 mjy, 
and 7 = -1.61 ± 0.01 and r = 5.21 ± 0.05ft _1 Mpc for 
Ssso/im ^ 1 mjy. For galaxies with 54501™ ^ 5 mjy, we find 
a slope similar to that of 850 /im galaxies, 7 — —1.62 ±0.01, 
and ro = 5.20T0.07 fr -1 Mpc. The values of ro obtained from 
the power law fit are thus close to the values obtained from 
the more general definition (Eq.[6]) in this case. However, we 
will use Eq. |S]to define ro unless stated otherwise. 

Our predictions for the clustering of the brighter SMGs 
at 850 fj,m are in reasonable agreem ent with the observa- 
tional estimate bv lBlain et al.l l|2004t ) who, using a sample of 
73 SMGs with iSsso/wn ^ 5 mjy and spectroscopic redshifts 
at z « 2 — 3, inferred a correlation length of 6.9T2.1 /i _1 Mpc 
using a pair-counting approach rather than a direct measure- 
ment of f(r) (note this measurement is discussed further in 
§3.7). 

We gain further insight into the clustering predicted by 
the model by plotting in Fig. [5] the mean number of sub- 
mm selected galaxies as a function of the halo mass, gener- 
ally refered to as th e halo occupation distribution or HOD 
dBerison et all l2000l : ICoorav fc Shethl |2002| : iBerlind et ail 
I2003T I. For completeness, we compute the median halo mass 
of our samples: for SMGs with Ssso/™ ^ 5 mjy and 
1 mjy, we find a median mass of 1.4 x 10 12 /i -1 Mq and 
9.3 x 10 11 /). -1 Mq, respectively; whereas for S^opm ^ 5 mjy 
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Figure 5. The predicted halo occupation distribution (HOD) of 
2 = 2 sub-mm selected galaxies. The dashed blue and orange 
lines show the HODs for galaxies brighter than 5 mjy at 450 and 
850 (im, while the green and red lines show the HODs for galaxies 
brighter than 1 mjy. The green dotted and dashed lines show 
the HOD of satellite and central galaxies with £450^™ ^ 1 mjy, 
respectively. 

and 1 mjy selected galaxies, we determine 9.3 x lO 11 ^" 1 M Q 
and 7.3 x lO 11 /! -1 M . Fig.[5]shows that both our 850 (im 
selected samples display a mean number of galaxies below 
unity over most of the range of halo masses. For example, in 
our model we find on average one Sssoum 5 mjy sub- 
mm galaxy for every ~ 100 dark matter haloes of mass 
~ 10 13 ft -1 M©. This illustrates the need to consider a large 
number of halo merger histories in order to make robust 
predictions, a point we made previously for luminous red 
galaxies (Almeida et al. 2008). We also predict that these 
massive haloes, with Mh a i 10 13 hT 1 Mq, will accommo- 
date more than one SMG with S^sopm 1 mjy (in our 
simulation, some haloes host more than 3 SMGs). If this 
was not the case, the two-point correlation function would 
tend to £ ~ —1 on s cales smaller than the typical size of 
the host haloes (see iBenson et all l2000h . For the case of 
5450/jm ^ 1 mjy we plot the contributions of central and 
satellite galaxies to the HOD separately. The central galaxy 
HOD is seen to flatten above a certain mass, while the satel- 
lite HOD continues to rise and becomes a power law. It is 
interesting to note in this case that the HOD for central 
galaxies does not reach unity for any halo mass; around 70% 
of haloes with masses 10 12 hT 1 Mq contain central galaxies 
brighter than 1 mjy at 450 /im. 

3.5 The redshift-space correlation function 

Galaxy surveys usually use redshift to infer the radial dis- 
tance to a galaxy, and the resulting measurement of clus- 
tering is said to be in redshift space. We therefore need 
to take into account the contribution of the peculiar ve- 
locity, induced by inhomogeneities in the galaxy's surround- 




-1 

log s (h _1 Mpc) 



Figure 6. The redshift space two-point correlation function of 
sub-mm selected galaxies. The correlation function for 450 fim 
and 850 (im selected galaxies with flux densities brighter than 5 
mjy, are represented by the dashed blue and orange lines, respec- 
tively, while solid green and red lines show galaxies with fluxes 
S450/jm ^ 1 rnjy and Sssopm ^ 1 mjy. The errorbars show the 
1<t Poisson errors derived from the number of pairs. For compar- 
ison we plot the S450jjm ^ 1 nijy real-space correlation function, 
using a dotted black line. 

ing density field, to the position of a galaxy inferred from 
its redshift. To model redshift space, we perturb the posi- 
tion of the galaxy along one of the cartesian axes, x, by 
the peculiar velocity of the galaxy along this axis, scaled by 
the appropriate value of the Hubble parameter. This corre- 
sponds to the redshift position as viewed by a distant ob- 
server. The redshift space correlation function for galaxies 
selected by their sub-mm flux is plotted in Fig. [6] The im- 
pact of including the peculiar velocities on the correlation 
function (redshift space distortions) depends on scale. On 
intermediate and large scales, the bulk motions of galaxies 
towards large scale structures generate an amplification of 
the amplitude of the correlation function (see the compari- 
son on large scales in I Jennings et"aT]|2010l ). On small scales, 
r < 3 h~ 1 Mpc, the peculiar motions of galaxies within struc- 
tures lead to a damping of the correlation function. In this 
case there is an apparent stretching of the structure in red- 
shift space, which dilutes the number of SMGs pairs. As a 
consequence, if we calculate the correlation length in redshift 
space, so, following the definition given by Eq. |6j i.e. so is 
given by £(so) = 1, we obtain: so = 6.4 ± 0.5 h~ Mpc for 
galaxies with Sssoum ^ 5 mjy, and so = 5.63 ±0.02 /i _1 Mpc 
for fainter galaxies Ssso^m ^ 1 mjy, slightly larger than the 
real-space values given in H3.4I 



3.6 The evolution of the correlation length 

Having computed the spatial correlation functions of SMGs 
at z = 2, it is natural to study the evolution with redshift of 
the clustering properties of these galaxies. In Fig. [7] we plot 
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Figure 7. Evolution of the comoving correlation length with 
redshift for galaxies selected with SssOfim 1? 1 rnjy (red line), 
•5850/jm ^ 5 mjy (dashed orange line), S 450^™ ^ 1 rnJy (solid 
green line), and S^so/jm ?S 5 mjy (dashed blue line). 



the evolution of the comoving real-space correlation length, 
ro, over the redshift interval 2 = — 4, for galaxies selected 
to be brighter than either 1 mjy or 5 mjy at 450 or 850 /im. 
Here we determine ro by finding the scale at which £(ro) = 1. 
The wavelengths quoted are in the observer's frame, which 
means that, for example, at z = 4 A = 850 /im corresponds 
to a rest-frame wavelength A rcst ~ 170 /im. Over the typical 
redshifts at which SMGs are found (i.e. around z = 2), Fig. [7] 
shows that the comoving correlation length is approximately 
constant. There is an increase in correlation length beyond 
z ~ 3 and a decrease at 2 < 1, but very few galaxies ap- 
pear in our samples at these redshifts. Fig. [7| indicates that 
the correlation length of galaxies selected at A = 450 /im 
is usually smaller than of 850 /im selected galaxies at the 
same flux limit. This is not surprising due to the fact that 
S45o A im/Ss50Mm ~ 3 (see the colour distribution in Fig. [T2")) , 
i.e. the same galaxy would appear about 3 times brighter at 
A = 450 /im than at A = 850 fim. This would then translate 
into a longer comoving correlation length (as we will see in 
Fig.©. 



3.7 Angular correlation function 

The simplest measure of clustering in photometric galaxy 
surveys is the angular two-point correlation function, w(8), 
which is a weighted projection of £(r) on the sky. The angu- 
lar correlation function is used to measure clustering when- 
ever redshift information is not available. In this case, the 
probability of finding two objects separated by angle 9 is 
similar to Eq. [2] 

6P{6) = ff [l + w(9)] SttiSQi, (8) 

where fj is the surface density of objects and 5£li is an ele- 
ment of solid angle. 



Figure 8. The two-point angular correlation function of sub-mm 
selected galaxies, computed directly from the spatial two-point 
correlation function. Galaxies selected at wavelengths of 450 fim 
and 850 /mi with fluxes brighter than 5 mjy are shown by dashed 
blue and orange lines, while the green and red lines show galaxies 
with S^soum ^ 1 m Jy an <l 5850jim 55 1 mjy respectively. The 
errorbars show the la Poisson errors derived from the number of 
pairs in each bin. These are computed for an area of 1.3 deg 2 at 
450 /xm and 35 deg 2 at 850 /im, to match upcoming surveys, as 
described in the text. 



The galaxy formation predicts the spatial two-point 
correlation function, £(r), and the redshift distribution of 
SMGs. From these predictions, it is straightforward to cal- 
culat e the angular correlation function using Limber's equa- 
tion l|Lmiberlll953r ). In a spatially flat universe, we have: 

f n °°du f n °° x 4 f 2 (x)ti(r,z)dx 
w(9)=2^ i2 V ; 3 ' , (9) 

[f °°x^(x)dx] 

where u is related to the comoving distance x by r 2 = 
u 2 + x 2 — 2xucos9. The selection function, ^(x), gives the 
probability that a sub-mm galaxy at a distance x is detected 
in the survey, and is defined by: 

pc<y -i poo 

K=/ x 2 <f(x)dx = — / N(z)dz, (10) 

JO Jo 

where H is the surface density of galaxies, f2 a the solid angle 
covered by the survey, and N(z) gives the number of sources 
within z and z + dz. 

In Fig. [8] we plot the angular correlation function for 
850 /im and 450 fim selected galaxies. We show predictions 
for galaxies brighter than 1 mjy and 5 mjy. We also show 
errorbars derived from the number of pairs in each bin in 
angular separation. These are computed for areas of 1.3deg 2 
and 35deg 2 at 450 fim and 850 /im respectively, chosen to 
match the survey areas planned with SCUBA-2. 

Fig. [8] shows that between 1 and 10 arcminutes , w(6) 
can be approximated by a power law: 
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where 7 is the same as in Eq. [7] The amplitude of clustering 
is characterised by 8q. This scale can vary from one sam- 
ple to another because of differences in the intrinsic cluster- 
ing and the redshift distribution of galaxies. Chance pro- 
jections of galaxy pairs, as would occur more often in a 
catalogue which spanned a broad redshift interval, lead to 
a dilution of the clustering signal, whi ch results in dp de- 
creasing with increasing survey depth (|Peebledll980h . The 
predicted angular correlation function steepens slightly be- 
low ~ larcmin. Fig. [8] shows that bright SMGs are more 
clustered than faint SMGs. Fitting Eq. [11] to our samples 
over the range of angular scales 1 < 8 < 10 arcmin, with a 
fixed 7 = —1.7 (see previous section), we find a clustering 
scale of (9o(<S850 M m ^ 5 mjy) = 0.028 ± 0.012 arcmin. For 
galaxies selected by A = 450 /im flux, we find 60(8450^ 5 
mjy) = 0.022 ± 0.007 arcmin. 

Observational estimates of the angular correlation func- 
tion of SMGs suggest that these galaxies are strongly clus- 
tered. However, a consensus on the clustering amplitude of 
SMGs is yet to be reached. For example, using a small sam- 
ple of SMGs observed with the SCUBA camera, IScott et al.l 
l|2002h found evidence of stro ng clustering on sc ales of 1 - 2 
arcmin. In a follow-up study. IScott et alT(|2006l ) re-reduced 
and combined different SCUBA surveys to measure the an- 
gular clustering and inferred a clustering scale 80 w 0.6 — 0.8 
arcmin, depending on the precise choice of flux limit and 
signal-to-noise cut-off used to construct maps of SMGs. The 
challenge of measuring the clustering signal with such small 
samples is apparent from the size of the integral constraint 
correction typically applied. The integral constraint takes 
into account the fact that the true mean density of SMGs is 
not known when estimating fluctuations in the galaxy dis- 
tribution. Instead, the density of the sample itself is used 
as a substitute for the true mean density. If the sample 
density is different from the true mean then the fluctua- 
tion level is misestimated, leading to a bias in w(8). This 
effect is important for small samples and is even more se- 
vere if the galaxies are also strongly clustered. Scott et al. 
applied a correction of a factor of » 2 to their measured clus- 
tering to account for the integral constraint. A subsequent 
measurement of clustering in a sample constructed with the 
LABOCA detector by Weiss et al. (2009) gives a correla- 
tion angle of 80 — 0.23 ± 0.12 arcmin, which is substantially 
weaker that the Scott et al. result. The Weiss et al. estimate 
is still much larger than our model prediction, but differs 
from it by only 1.8<r. Blain et al. (2004) made an indirect 
measurement of the spatial correlation function using radial 
pair counts of SMGs with spectroscopic redshifts, and ob- 
tained ro = 6.9 ± 2.1ft _1 Mpc. For spatial clustering that is 
fixed in comoving coordinates, and using the observed red- 
shift distribution of SMGs (Chapman et al. 2004), the Blain 
et al. result implies 80 — 0.04 ± 0.01 arcmin, which is com- 
parable to our model prediction However, the method used 
by Blain et al. does not take into account fluctuations in the 
sample density, and should perhaps be viewed as providing 
a lower limit on the clustering strength. 
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Figure 9. The dependence of the comoving correlation length, 
ro, on sub-mm flux for galaxies at z = 2. The predictions for 
SMGs selected by their emission at 850 fim are shown by the red 
line, and those for 450 (im selected sub-mm galaxies are plotted in 
blue. The flux plotted on the x-axis is correspondingly at either 
850 or 450 fj,m. 



4 THE DEPENDENCE OF CLUSTERING 
STRENGTH ON GALAXY PROPERTIES 

Here we consider the dependence of clustering strength on 
galaxy properties. We present the predictions for the two- 
point correlation function as a function of sub-mm flux, halo 
and stellar mass, S45o Mm /S850Mm colour and the quiescent or 
starburst nature of galaxies at z — 2. 



4.1 Dependence of clustering on sub-mm flux 

In Fig. [9] we plot the comoving correlation length, ro (com- 
puted using Eq. 6), as a function of sub-mm flux, S„, for 
galaxies at 2 = 2. The plot shows that brighter galaxies 
have larger correlation lengths (see also Figs.|4]and [8]). How- 
ever, the dependence of clustering strength on luminosity is 
fairly weak, with a change of 50% in correlation length on 
changing flux by a factor of a hundred. This behaviour can 
be understood as a consequence of brighter galaxies being 
found predominately in more massive haloes, for which the 
bias is greater than unity and increases strongly with mass 
(Angulo et al. 2009). As already pointed out in relation to 
Fig. [4] for the same flux limit, 450 selected galaxies are 
less clustered than their 850 /im counterparts, with correla- 
tion lengths that are typically smaller by Aro ~ 0.4 /i _1 Mpc. 
This difference remains roughly constant throughout the 
range of sub-mm fluxes explored (in Fig. [T3] we will see 
that this is a consequence of the fact that the median 
S450Mm/5 , 850Mm colour is approximately w 3). 
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Figure 10. The dependence of the comoving correlation length, 
ro, on host halo mass, Afhaloj for galaxies at z = 2 (for differential 
bins in halo mass). The red and green solid lines display the rela- 
tion predicted for galaxies with Ssso^m ^ 1 m Jy an d S4so^m 1 
mjy respectively. Brighter sub-mm galaxies, with fluxes above 
5 mjy, are shown by the dashed orange (850 fim) and dashed 
blue (450 fim) lines. The dashed black line shows the dependence 
of correlation length on halo mass for all haloes (regardless of 
whether they host an SMG or not). 
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Figure 11. The dependence of the comoving correlation length, 
ro, on total stellar mass, A/,, for galaxies at z = 2. The red 
and green solid lines display the predictions for galaxies with 
SssO/jm 1 mJ y and S450/*m 3= 1 mjy, respectively. Brighter 
sub-mm galaxies, with fluxes brighter than 5 mjy, are shown 
by the dashed orange (850 fim) and dashed blue (450 fim) lines. 
The relation between ro and stellar mass for all the galaxies in 
the simulation, regardless of their sub-mm flux, is shown by the 
dashed-black line. 



4.2 Dependence of clustering strength on halo 
and stellar mass 

Fig. [10] shows the dependence of ro on the mass of the dark 
matter halo which hosts the SMG, for discrete bins in halo 
mass (the bin width is adjusted to contain a representative 
sample of galaxies). Remarkably, there is little difference be- 
tween the clustering strength predicted for bright and faint 
samples of SMGs within a given bin of host halo mass. More- 
over, the SMG clustering is essentially the same as that of 
all of the dark matter haloes in the mass bin, regardless of 
whether or not they contain an SMG. This means that the 
clustering strength of SMGs is driven purely by the mass of 
the host halo, and shows little or no dependence on a second 
property of the halo, such as its formation time or spin. One 
might have expected to see a dependence of the clustering 
strength at a given halo mass on SMG flux if this selection 
favoured host haloes which had, for example, formed more 
rece ntly than the overall population at a given mass (see, 
e.g. iGao et al.|[2005l : iPercival et a"i1l2003h . 

We now consider the dependence of the correlation 
length ro on stellar mass, M*, which we plot in Fig. 1111 Fol- 
lowing earlier plots, we display the relation for SMGs with 
fluxes brighter than 1 mjy at 450/xm and 850pim, and for 
brighter samples with fluxes ^ 5 mjy at both wavelengths. 
At z = 2, the median stellar masses of the model samples 
are M* = 9.9 x 10 9 /i _1 Mq for galaxies with 5*850 > lmJy, 
and Aft = 9.7 x W 9 h~ 1 Mp for those with S450 > lmJy (see 
also lGonzalez et alj|2010bh . Comparisons with observational 
estimates of stellar masses are of limited use, as the results 
depend critically on the assumption made about the form 
of the stellar IMF (see La cey et al. 2010a for a n expanded 
discussion). For example, lHainline et al.l (|2010l ) estimated 
the stellar masses of ~ 70 SMGs with Sssotim ~ 5 mjy and 
found a median stellar mass of 7 x 10 10 h- 1 M e , when using 
a Kroupa (2001) IMF. At the same flux limit, our model 
predicts a median stellar mass of 2 x 1Q 10 /i- 1 M q . 

There is a tight, monotonic correlation between the co- 
moving correlation length, ro, and stellar mass, M* (Fig. llip . 
Massive galaxies are more strongly clustered than galaxies 
with lower stellar masses, implying that more massive galax- 
ies tend to be hosted by more massive haloes. Interestingly, 
all of our samples show a similar correlation length for a 
given total stellar mass, i.e. the relation between correlation 
length and total stellar mass does not depend on the sub- 
mm luminosity of the galaxy. We find that the correlation 
length ro varies roughly linearly with log M* in the range 
log Af, ~ [9.5, 10.5]/i _1 Mq. Another important aspect is the 
fact that SMGs have, typically, weaker clustering than other 
galaxies of the same stellar mass. This is particularly evident 
for stellar masses below 10 10 /i -1 Mq. The phenomenon is es- 
sentially due to the combination of two factors. Firstly, there 
is a difference in the host halo masses for satellite and cen- 
tral galaxies of the same stellar mass: at a given stellar mass, 
satellite galaxies inhabit haloes which are roughly ten times 
more massive than central galaxies. Secondly, SMGs galax- 
ies with stellar masses lower than lO lo /i~ 1 M0 are mainly 
central galaxies, whereas the fraction of satellite galaxies for 
the full galaxy population is approximately 30 per cent. Sub- 
millimetre galaxies with stellar masses 
display a similar fraction of satellite galaxies to that in the 
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Figure 12. Distribution of the 8450^111 /SssOum ^ ux ratio or 
colour for SMGs at 2 = 2, as a function of 850 (im flux density. 
The solid red and dashed orange histograms show the distribution 
for quiescent and starburst galaxies, respectively, while the black 
line displays the combined distribution. Quiescent galaxies only 
make a significant contribution at the faintest fluxes (top left). 
The distributions are normalized to integrate to unity. 

full sample. Currently, there are no observational determi- 
nations of this relation at the redshifts of interest for SMGs. 



4.3 Dependence of clustering on sub-mm colour 

Before studying the dependence of the correlation length, r , 
on sub-mm colour, as given by the S450(jm/Ss50Mm flux ratio, 
it is informative to first plot the colour distribution itself, 
which we show in Fig. 1121 Here, we separate galaxies into 
bins of S850 M m flux and discriminate between quiescent and 
starburst galaxies. In our model, the S450/im/S850jxm colour 
distribution of A = 850 /im selected galaxies peaks around 
« 3.2. Furthermore, we find that the mode and shape of the 
colour distribution does not change significantly with Ssso^m 
flux, i.e. the colour-luminosity relation of sub-mm galaxies 
with Sgso/jm ^ 1 rnjy is roughly constant. Also, there is 
very little difference between the colour distribution of star- 
burst and quiescent galaxies (note, however, that, at lower 
fluxes, the fraction of sub-mm galaxies which are forming 
stars quiescently increases). 

The observed distribution of 450/im/850/im colours 
of SMGs has not yet been accurately measured, but the 
350/im/850At m colours of SMG s in th e same model were in- 
vestigated by S winbank et al l <|2008l ). who found the pre- 
dicted colours to be in agreement with observed values. 
They also found that the median 350/^m/850/im colour of 
model SMGs could be fit by a modified blackbody spectrum 
L v oc Bv(T) v 13 , with /? = 1.5 and an effective dust tempera- 
ture of T — 32 K. This modified blackbody implies a colour 
5 , 450Aim/Sg5o A im ~ 3.5 for SMGs at z = 2, which agrees well 
with the distribution plotted in Fig. 1131 

In Fig. [13] we plot the correlation length, ro, as a func- 




log S x (mJy) 



Figure 13. Dependence of the correlation length, ro, on 
S4B0fim / SsBOfim colour for galaxies at z = 2. SMGs selected by 
their emission at A = 850 /an are represented by the red line, while 
the predictions for galaxies selected at A = 450 fira are plotted in 
blue. Solid and dashed lines show the correlation length for galax- 
ies split by sub-mm colour, S45o M m/'S l 85o t im> as indicated by the 
key. 

tion sub-mm flux for two samples split at a sub-mm colour of 
S450fim/S850fim = 3.2. The dependence of ro on luminosity is 
similar to that shown in Fig. [9] The figure also hints that red- 
der galaxies, i.e. those with S45o^m/S850,jm ^ 3.2, are more 
clustered than bluer galaxies with S450Mm/S850/im ^ 3.2 (for 
both A = 850 /xm and A = 450 [im selected galaxies). For 
example, for Ssso/jm ~ 3 mJy galaxies, ro can differ by a 
factor of ~ 1.6 between the red and blue samples. 

Having plotted the colour distributions of starburst and 
quiescent galaxies and the relation between ro and sub-mm 
colour, it is useful to study the clustering of quiescent and 
starburst galaxies separately. As mentioned in connection 
with Fig. 1121 quiescent galaxies only make a significant con- 
tribution at lower fluxes, SsBo^m ^ 2 mJy. In Fig. 1141 we plot 
the two-point correlation function in real-space of sub-mm 
galaxies selected by Sssom™ ^ 1 m Jy at z = 2. We also plot 
the correlation functions of quiescent and starburst galax- 
ies separately. Fig. 1131 reveals that, in our model, quiescent 
galaxies are more clustered than burst galaxies of similar 
sub-mm luminosities, on all scales. This is mostly due to 
the fact that quiescent SMGs are more massive than burst 
SMGs, due to both the top-heavy IMF and shorter star for- 
mation timescale in bursts. We compute ro = 5.2 h~ x Mpc 
for burst galaxies with Sssoum 1 rnjy, and ro = 7.2 h~ l 
Mpc for quiescent galaxies brighter than the same flux limit. 



5 DISCUSSION AND CONCLUSIONS 

In this paper we hav e applied the technique introduced by 
lAlmeida et"afl (|2010h . based on artificial neural networks, 
to predict the spatial and angular clustering of sub-mm se- 
lected galaxies (SMGs) in a ACDM universe. The ANN al- 
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Figure 14. The real-space two-point correlation function of qui- 
escent (red line) and starburst (orange) sub-mm galaxies, selected 
by Sgsopm i? 1 rnJy at z = 2. The dashed black line shows the 
correlation function for all galaxies with fluxes brighter than 1 
mjy. 



lows us to rapidly mimic the predictions of a hybrid code, 
made up of the GALFORM semi-analytical model of galaxy 
formation, which predicts the full star formation and merger 
histories of galaxies, and the GRASIL code, which computes 
the spectral energy distributions of galaxies self-consistently. 
This makes it possible for large N-body simulations of the 
hierarchical growth of structure in the dark matter to be 
populated with galaxies with full SED coverage from the UV 
through to the radio. The ANN technqiue requires as input 
a small set of galaxy properties predicted by GALFORM. We 
use the algorithm to populate the Millennium Simulation 
with A = 850 /im and A = 450 /im selected galaxies, with 
3\ ^ 1 mjy. The accuracy of the artificial neural network 
is notable: for SMGs at z = 2, we are able to reproduce the 
luminosity of ~ 95 per cent of galaxies to within 10 per cent 
of the true values. The errors introduced by this technique 
in the determination of clustering properties are negligible. 

We have presented predictions for the two-point spa- 
tial and angular correlation functions, for samples of galax- 
ies selected at 850 /im and 450 /mi. At z = 2, we predict a 
comoving correlation length of ro = 5.6 ± 0.9 /i -1 Mpc for 
galaxies with Sgso^m 5 mjy and ro = 5.38 ± 0.02 /i _1 Mpc 
for fainter galaxies, brighter than 1 mjy. The former value is 
in good agreement w ith the indirect observational estimate 
bv lBlain etakl (|2004T ) who found r = 6.9 ± 2.1 /i _1 Mpc for 



a small sample of SMGs with SssOfj 



5 mjy. Galaxies se- 



lected at A = 450 /im are less clustered than those identified 
at 850 /im galaxies, to the same flux limit. 

Not surprisingly, the correlation length of sub-mm se- 
lected galaxies evolves with redshift: galaxies selected at 
higher redshifts display larger comoving correlation lengths, 
i.e. with increasing redshift they trace more overdense (and 
consequently, rarer) regions of the Universe. This behaviour 
is similar to that expected for the growth of structures in a 



ACDM universe, and further supports the idea that sub-mm 
galaxies are biased tracers of the underlying mass. 

We have also studied the dependence of clustering 
on some properties of our sub-mm galaxies. Generally, we 
found a strong dependence of clustering on luminosity, with 
brighter galaxies displaying higher correlation lengths. We 
predict a tight correlation between halo mass and clustering 
scale, with sub-mm galaxies hosted by massive haloes being 
more clustered than those hosted by less massive ones. We 
find that the relation found for SMGs is not significantly 
different from the clustering of all haloes of the same mass. 
A similar behaviour is observed for the correlation between 
clustering and total stellar mass. We find that more mas- 
sive sub-mm galaxies are more clustered than less massive 
one (in fact, massive galaxies are predominately found in 
more massive haloes - haloes with higher effective bias): in 
our model, ro changes by a factor of ~ 2 within the range 
log Alt ~ [9.5, 10.5] h~ Mq. The dependence of clustering on 
sub-mm colour, SAbo^m/Ssso^m, is not so clear. We find a 
weak positive correlation between ro and colour, with redder 
galaxies being more strongly clustered. Finally, we predict 
that for galaxies selected to have Ssso^m ^ 1 mjy, quiescent 
galaxies are more clustered than those undergoing a burst 
of star formation. 

It is important to note that current observational mea- 
surements of the clustering properties of sub-mm selected 
galaxies are limited due to poor statistics. At z ~ 2, estima- 
tions of the correlation function currently rely on samples 
with fewer than ~ 100 galaxies. However, this picture will 
soon be improved. A more detailed comparison will soon 
be possible with the forthcoming instruments, such as the 
SCUBA-2, which will allow a much deeper and wider sur- 
vey of the submillimetre population up to redshift 4 than 
current instruments. A more accurate determination of the 
clustering will, by then, impose serious constrains on the 
theories of galaxy formation and evolution. 
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